package fem2;

import inf.math.LinAlg;
import math2.IsoparametricFaceGeometry;
import math2.LinAlgUtilities;
import math2.MathUtilities;
import fem2.enu.MeshPartType;

/**
 * 
 * TODO
 * 
 * @author hbui
 * 
 */
public class Face extends MeshPart {

	private Node[] nodes;

	/**
	 * constructor for isoparametric face
	 * 
	 * @param nodes
	 */
	public Face(Node... nodes) {
		this.nodes = nodes;
		if (nodes.length == 3) {
			setType(MeshPartType.TRI3);
		} else if (nodes.length == 4) {
			setType(MeshPartType.QUAD4);
		} else if (nodes.length == 6) {
			setType(MeshPartType.TRI6);
		} else if (nodes.length == 8) {
			setType(MeshPartType.QUAD8);
		} else {
			throw new Error("this geometry is not supported");
		}
		setBasis(BasisManager.getInstance().getBasis(getType()));
		setGeometry(IsoparametricFaceGeometry.getInstance());
		/*
		 * chec the positiveness of jacobian to ensure correct mapping
		 */
		if (getDimension() == 2) {
			boolean posJ = checkJacobian();
			if (!posJ) {
				String ids = "";
				int nn = nodes.length;
				for (int i = 0; i < nn; i++) {
					ids += nodes[i].getId() + " ";
				}
				// throw new Error(
				// "negative jacobian is not allowed, check face with nodes "
				// + ids);
				System.out.println("warning: negative jacobian is detected at face " + ids);
				System.out.println("modified jacobian by swapping node indices");
				this.nodes = new Node[nn];
				for (int i = 0; i < nn; i++) {
					this.nodes[i] = nodes[nn - 1 - i];
				}
				System.out.print("recheck jacobian...");
				posJ = checkJacobian();
				if (!posJ) {
					throw new Error("still negative. The connectivity is invalid.");
				} else {
					System.out.print("passed\n");
				}
			}
		}
	}

	@Override
	public Node[] getNodes() {
		return nodes;
	}

	@Override
	public double[][] gGradientsAt(double... xi) {
		if (getDimension() == 2) {
			double[][] iJac = this.jacXAt(xi);
			LinAlg.invert2x2(iJac);
			int n = this.getBasis().getDimension();
			double[][] gGradient = new double[n][2];
			for (int i = 0; i < n; i++) {
				double[] tmp = this.getBasis().getBasisFunction(i).gradientAt(xi);
				double[] dNi = new double[2];
				LinAlg.multiply(false, 2, 2, 1.0, iJac, tmp, 0.0, dNi);
				gGradient[i][0] = dNi[0];
				gGradient[i][1] = dNi[1];
			}
			return gGradient;
		} else if (getDimension() == 3) {
			/*
			 * TODO
			 */
			return null;
		} else {
			throw new Error("invalid dimension");
		}
	}

	/**
	 * compute area in generic way. The difference with compuyeVolume is that
	 * computeArea can compute area in 3d
	 * 
	 * @return
	 */
	public double computeArea() {
		// return computeVolume();
		IntegrationManager integrator = IntegrationManager.getInstance();
		double[] xi = integrator.getIntegrationPoint(getType(), 1, 0);
		double w = integrator.getIntegrationWeight(getType(), 1, 0);
		double[][] jac = jacXAt(xi);
		return w * Math.sqrt(LinAlgUtilities.determinantJJt(jac));
	}

	@Override
	public boolean checkInside(double... p) {
		assert p.length == getDimension() : "incompatible dimension";

		double area = 0.0;
		for (int i = 0; i < nodes.length; i++) {
			Node n1 = nodes[i];
			Node n2;
			if (i == nodes.length - 1) {
				n2 = nodes[0];
			} else {
				n2 = nodes[i + 1];
			}
			area += MathUtilities.heron(p, n1.getPosition(), n2.getPosition());
		}

		return Math.abs(area - computeArea()) < MathUtilities.tol6;
	}

}
